Auxiliary functions
get_deseq_data <- function(target_tissue, meta_dt, data_dir, tx2knownGene) {
samples_dir <- file.path(
data_dir,
meta_dt %>%
pull(name)
)
cdna_dir <- file.path(samples_dir)
cdna_files <- file.path(cdna_dir, "abundance.h5")
names(cdna_files) <-
meta_dt %>%
pull(sample)
txi_scaledTPM <- tximport(cdna_files,
type = "kallisto",
tx2gene = tx2knownGene,
countsFromAbundance = c("scaledTPM")
)
txi <- tximport(cdna_files,
type = "kallisto",
tx2gene = tx2knownGene,
countsFromAbundance = c("no")
)
scaled_tpm <-
txi_scaledTPM$counts %>%
as_tibble(rownames = "gene_id")
gene_counts <-
txi$counts %>%
as_tibble(rownames = "gene_id")
write_tsv(scaled_tpm, paste0(target_tissue, "_engraftment_scaled_tpm.tsv"))
write_tsv(gene_counts, paste0(target_tissue, "_engraftment_counts.tsv"))
design_matrix <-
model.matrix(
~ time * age * batch,
meta_dt %>%
column_to_rownames("sample")
)
colnames(design_matrix) <- c(
"intercept",
"engraf",
"age",
"intercept_batch",
"niche",
"engraf_batch",
"age_batch",
"niche_batch"
)
design_matrix <- design_matrix[, c(1, 2, 3, 5, 4, 6, 7, 8)]
dds <-
DESeqDataSetFromTximport(
txi,
meta_dt %>%
column_to_rownames("sample"),
design = design_matrix
)
return(dds)
}
# We use apeglm shrinkage to preserve the size of large LFC and compute s-values (the estimated rate of false sign). Shrinkage is useful for ranking and visualization, without the need for arbitrary filters on low count genes. https://doi.org/10.1093/bioinformatics/bty895
# The local false sign rate FSR is defined as the posterior probability that the posterior mode (MAP) of beta (log2FC) is of a biologically significant effect size (i.e. abs(beta) > lfcThreshold )
get_lfc <- function(dds, coef_name, lfc_threshold, svalue_threshold) {
lfc <-
lfcShrink(dds,
coef = coef_name,
type = "apeglm",
svalue = T,
lfcThreshold = lfc_threshold
)
plotMA(lfc, ylim = c(-22, 8))
abline(h = c(-1, 1), col = "dodgerblue", lwd = 2)
lfc <-
lfc %>%
as_tibble() %>%
dplyr::select(-baseMean) %>%
mutate(
keep = ((abs(log2FoldChange) > lfc_threshold) & (svalue < svalue_threshold))
) %>%
set_names(c(
paste0(
c("lfc_", "lfc_se_", "svalue_", "keep_"),
coef_name
)
))
return(lfc)
}
get_lfc_estimates <- function(dds, gene_metadata, coef_names, lfc_threshold, svalue_threshold, batch_threshold) {
lfcs <-
map_dfc(coef_names,
get_lfc,
dds = dds,
lfc_threshold = lfc_threshold,
svalue_threshold=svalue_threshold
)
coef_names_batch <- c(paste0(c("intercept", coef_names), "_batch"), paste0("SE_", c("intercept", coef_names), "_batch"))
lfcs <-
bind_cols(
mcols(dds)[, c("baseMean", "baseVar", "dispOutlier" , "maxCooks", "intercept", "SE_intercept", c(coef_names, coef_names_batch))] %>%
as_tibble(rownames = "gene_id") %>%
dplyr::select(gene_id, everything()),
lfcs
)
lfcs <-
lfcs %>%
mutate(
total_batch= abs(intercept_batch) + abs(niche_batch) + abs(age_batch) +abs(engraf_batch),
keep_niche = (keep_niche & abs(niche_batch) < lfc_threshold & total_batch < batch_threshold ),
keep_age = (keep_age & abs(age_batch) < lfc_threshold & total_batch < batch_threshold ),
keep_engraf = (keep_engraf & abs(engraf_batch) < lfc_threshold & total_batch < batch_threshold ),
EAN = paste0(
as.integer(keep_engraf),
as.integer(keep_age),
as.integer(keep_niche)
),
)
lfcs <-
lfcs %>%
dplyr::select(
gene_id, EAN, baseMean,
intercept, paste0("lfc_", coef_names),
paste0("svalue_", coef_names), everything()
)
lfcs <-
lfcs %>%
left_join(., gene_metadata, by = "gene_id") %>%
dplyr::select(gene_name, everything()) %>%
mutate(
gene_name =
scater::uniquifyFeatureNames(
ID = gene_id,
names = gene_name
)
)
return(lfcs)
}
plot_volcano <- function(data, xvar, yvar, title, lfc_threshold, svalue_thresh, figures_dir) {
gv <-
EnhancedVolcano(data,
lab = data %>% pull(gene_name),
xlab = bquote(~ italic(Moderated) ~ Log[2] ~ FC),
ylab = bquote(~ -Log[10] ~ italic(svalue)),
x = xvar,
y = yvar,
col = c("grey30", "forestgreen", "red2", "royalblue"),
pCutoff = svalue_thresh,
FCcutoff = lfc_threshold,
# legend = c("NS", "Log2FC", "svalue", "svalue & Log2FC"),
legendLabels = c(
"NS",
expression(Log[2] ~ FC),
"svalue",
expression(s - value ~ and ~ log[2] ~ FC)
)
)
ggsave(file.path(figures_dir, paste0(title, "_volcano.pdf")), gv)
return(gv)
}
plot_heatmap <- function(target_genes, title, meta_col, data, show_genes = F, figures_dir) {
df <-
left_join(target_genes %>%
dplyr::select(gene_id, gene_name),
data,
by = "gene_id"
) %>%
dplyr::select(-gene_id) %>%
column_to_rownames("gene_name") %>%
as.matrix()
gh <-
pheatmap(df,
color = colorRampPalette(rev(RColorBrewer::brewer.pal(n = 11, name = "Spectral")))(100),
cluster_rows = T,
show_rownames = show_genes,
treeheight_row = 0,
fontsize_col = 4,
fontsize_row = 6,
angle_col = "45",
cluster_cols = FALSE,
annotation_col = meta_col
)
ggsave(file.path(figures_dir, paste0(title, "_heatmap_expression.pdf")), gh)
return(gh)
}
plot_pca <- function(vsd, title, figures_dir) {
pcaData <-
plotPCA(vsd,
intgroup = c("time", "age", "tissue", "batch", "donor"),
returnData = TRUE
)
percentVar <- round(100 * attr(pcaData, "percentVar"))
gg_plot <-
ggplot(
pcaData %>%
mutate(condition = paste(age, time, sep = "_")),
aes(x = PC1, y = PC2, color = time, shape = age, label = name)
) +
geom_point(size = 4) +
geom_line(aes(group = donor), colour = "grey") +
# coord_fixed() +
geom_text(check_overlap = T, angle = 0, size = 3, vjust = 3, nudge_y = 1.5) +
xlab(paste0("PC1: ", percentVar[1], "% variance")) +
ylab(paste0("PC2: ", percentVar[2], "% variance")) +
coord_fixed() +
ggtitle("PCA with VST data")
ggsave(file.path(figures_dir, paste0(title, "_pca_bulk_rnaseq.pdf")), gg_plot)
gg_plot
}
plot_gene <- function(gene, data, figures_dir) {
df_plot <-
data %>%
dplyr::filter(gene_id == gene)
gene_name <- df_plot$gene_name[1]
p <- ggplot(
df_plot,
aes(x = time, y = log2_scaled_tpm)
) +
geom_point(aes(color = donor)) +
geom_line(aes(group = donor)) +
labs(title = gene_name) +
facet_grid(batch ~ age)
ggsave(file.path(figures_dir, paste0(gene_name, "_expression_log2_scaled_tpm.pdf")), p)
return(p)
}
Get tx id versions
library(AnnotationHub)
Loading required package: BiocFileCache
Loading required package: dbplyr
Attaching package: 'dbplyr'
The following objects are masked from 'package:dplyr':
ident, sql
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Attaching package: 'AnnotationHub'
The following object is masked from 'package:Biobase':
cache
ah <- AnnotationHub()
DEPRECATION: As of AnnotationHub (>2.23.2), default caching location has changed.
Problematic cache: /home/rick/.cache/AnnotationHub
See https://bioconductor.org/packages/devel/bioc/vignettes/AnnotationHub/inst/doc/TroubleshootingTheCache.html#default-caching-location-update
snapshotDate(): 2021-05-18
meta_dt <-
read_tsv(file.path(metadata_filepath)) %>%
# arrange(tissue, age, time) %>%
mutate(
condition = paste(age, time, sep = "_"),
age = factor(age, levels = c("young", "aged")),
time = factor(time, levels = c("T0", "T21")),
batch = factor(batch, levels = c(1, 2))
)
Rows: 23 Columns: 8── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (7): name, sample, time, tissue, age, donor, type
dbl (1): batch
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
meta_dt
gene_metadata <-
read_csv(gene_metadata_filepath) %>%
dplyr::select(-gene_id_version)
Rows: 55487 Columns: 8── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (6): gene_name, gene_id, seqnames, gene_biotype, gene_id_version, description
dbl (2): gene_len, perc_gene_gc
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
gene_metadata <-
left_join(tx_meta %>% as_tibble() %>%
group_by(gene_id) %>%
dplyr::count(gene_name) %>%
dplyr::select(-n),
gene_metadata %>%
dplyr::select(-gene_name),
by = "gene_id"
) %>%
mutate(
gene_name =
scater::uniquifyFeatureNames(
ID = gene_id,
names = gene_name
)
)
Load data
dds_leg <- get_deseq_data("leg", meta_dt, data_dir, tx2knownGene)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
summarizing abundance
summarizing counts
summarizing length
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
summarizing abundance
summarizing counts
summarizing length
using counts and average transcript lengths from tximport
dds_leg
class: DESeqDataSet
dim: 36558 23
metadata(1): version
assays(2): counts avgTxLength
rownames(36558): ENSMUSG00000000001 ENSMUSG00000000003 ... ENSMUSG00000118487 ENSMUSG00000118488
rowData names(0):
colnames(23): yng_t0_01 yng_t0_02 ... aged_t21_05 aged_t21_06
colData names(8): name time ... batch condition
summary_dt_leg <-
colSums(assay(dds_leg)) %>%
tibble::enframe(name = "sample", value = "counts")
summary_dt_leg
scaled tpm
scaled_tpm_leg <- read_tsv("leg_engraftment_scaled_tpm.tsv")
Rows: 36558 Columns: 24── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (1): gene_id
dbl (23): yng_t0_01, yng_t0_02, yng_t0_03, yng_t0_04, yng_t0_05, yng_t0_06, yng_t21_01, yng_t21_02, yng_t21_04, yng_t21_05, yng_t21_06, aged_t0_01, aged_t0_02, aged_t0_03, aged_t0_04, aged_t0_05, aged_t0_06, aged_t21_01, aged_t21_02, a...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
counts_leg <- read_tsv("leg_engraftment_counts.tsv")
Rows: 36558 Columns: 24── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (1): gene_id
dbl (23): yng_t0_01, yng_t0_02, yng_t0_03, yng_t0_04, yng_t0_05, yng_t0_06, yng_t21_01, yng_t21_02, yng_t21_04, yng_t21_05, yng_t21_06, aged_t0_01, aged_t0_02, aged_t0_03, aged_t0_04, aged_t0_05, aged_t0_06, aged_t21_01, aged_t21_02, a...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
scaled_tpm_leg_long <-
scaled_tpm_leg %>%
pivot_longer(cols = yng_t0_01:aged_t21_06, names_to = "sample", values_to = "scaled_tpm") %>%
left_join(., meta_dt %>% dplyr::select(sample, time, age, donor, batch, condition), by = "sample")
scaled_tpm_leg_long <-
left_join(scaled_tpm_leg_long,
gene_metadata %>%
dplyr::select(gene_name, gene_id),
by = "gene_id"
) %>%
mutate(log2_scaled_tpm = log1p(scaled_tpm))
summary_tpm_leg_long <-
scaled_tpm_leg_long %>%
group_by(gene_id, condition) %>%
summarize(
mean_tpm = mean(log2_scaled_tpm),
sd_tpm = sd(log2_scaled_tpm),
frac_zeros = mean(scaled_tpm == 0)
) %>%
mutate(cv = sd_tpm / mean_tpm) %>%
group_by(gene_id) %>%
mutate(
low_frac_zeros = sum(frac_zeros <= 0.2),
sd_group_tpm = sd(mean_tpm),
mean_group_tpm = mean(mean_tpm),
min_tpm = min(mean_tpm),
max_tpm = max(mean_tpm)
) %>%
dplyr::filter(low_frac_zeros > 0) %>%
# dplyr::filter(min_tpm >0.1) %>%
dplyr::filter(max_tpm > log(30))
`summarise()` has grouped output by 'gene_id'. You can override using the `.groups` argument.
summary_tpm_leg_long
ggplot(summary_tpm_leg_long %>%
group_by(gene_id) %>%
dplyr::slice(1)) +
geom_point(aes(x = min_tpm, y = max_tpm - min_tpm), size = 0.5)

target_genes <-
summary_tpm_leg_long %>%
pull(gene_id) %>%
unique()
keep <- rownames(counts(dds_leg)) %in% target_genes
dds_leg <- dds_leg[keep, ]
Generalized Principal Component Analysis
set.seed(11678)
gpca <- glmpca(scaled_tpm_leg %>%
dplyr::filter(gene_id %in% target_genes) %>%
column_to_rownames("gene_id"),
L = 2,
optimizer = c("fisher"),
fam = c("poi"),
X = model.matrix(~ -1 + batch, meta_dt %>% dplyr::select(batch))
)
gpca.dat <- gpca$factors
gpca.dat <- cbind(gpca.dat, meta_dt[, c("age", "time", "batch", "donor")])
gpca.dat
gg_pca <-
ggplot(
gpca.dat %>% as_tibble(rownames = "sample") %>%
mutate(condition = paste(age, time, sep = "_")),
aes(x = dim1, y = dim2, color = condition, shape = batch, label = sample)
) +
geom_point(size = 4) +
geom_line(aes(group = donor), colour = "grey") +
# coord_fixed() +
geom_text(check_overlap = T, angle = 0, size = 3, vjust = 2, nudge_y = 1.5) +
# facet_grid(~tissue) +
ggtitle("glmpca - Generalized PCA")
gg_pca

colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
poisd <- PoissonDistance(t(counts(dds_leg)))
samplePoisDistMatrix <- as.matrix( poisd$dd )
rownames(samplePoisDistMatrix) <- colnames(counts(dds_leg))
colnames(samplePoisDistMatrix) <- colnames(counts(dds_leg))
pheatmap(samplePoisDistMatrix,
clustering_distance_rows = poisd$dd,
clustering_distance_cols = poisd$dd,
col = colors)

vsd <- vst(dds_leg, blind = T)
using 'avgTxLength' from assays(dds), correcting for library size
pcaData <- plotPCA(vsd, intgroup = c("condition", "batch"), returnData = TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(x = PC1, y = PC2, color = condition, shape = batch)) +
geom_point(size =3) +
xlab(paste0("PC1: ", percentVar[1], "% variance")) +
ylab(paste0("PC2: ", percentVar[2], "% variance")) +
coord_fixed() +
ggtitle("PCA with VST data")

design_matrix <-
model.matrix(
~1,
meta_dt %>%
column_to_rownames("sample")
)
colnames(design_matrix) <- c("intercept")
dds_leg <- DESeq(dds_leg, test = "LRT", reduced = design_matrix)
using supplied model matrix
estimating size factors
using 'avgTxLength' from assays(dds), correcting for library size
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
plotDispEsts(dds_leg)

res <- results(dds_leg)
W <- res$stat
maxCooks <- apply(assays(dds_leg)[["cooks"]],1,max)
idx <- !is.na(W)
plot(rank(W[idx]), maxCooks[idx], xlab="rank of Wald statistic",
ylab="maximum Cook's distance per gene",
ylim=c(0,5), cex=.4, col=rgb(0,0,0,.3))
m <- ncol(dds_leg)
p <- 3
abline(h=qf(.99, p, m - p))

plot(res$baseMean+1, -log10(res$pvalue),
log="x", xlab="mean of normalized counts",
ylab=expression(-log[10](pvalue)),
ylim=c(0,30),
cex=.4, col=rgb(0,0,0,.3))

lfc_threshold <- 1
svalue_thresh <- 0.05
batch_threshold <- 4
coef_names <- c("engraf", "age", "niche")
lfcs <-
get_lfc_estimates(dds_leg,
gene_metadata,
coef_names,
lfc_threshold,
svalue_thresh,
batch_threshold)
using 'apeglm' for LFC shrinkage. If used in published research, please cite:
Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
sequence count data: removing the noise and preserving large differences.
Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)
thresholding s-values on alpha=0.005 to color points

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
sequence count data: removing the noise and preserving large differences.
Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)
thresholding s-values on alpha=0.005 to color points


lfcs
Plot
p1 <- ggplot(
lfcs,
aes(x = intercept_batch, y = intercept)
) +
geom_point(size = .7, alpha = 0.3) +
geom_vline(xintercept = c(-lfc_threshold, lfc_threshold)) +
geom_hline(yintercept = c(-1, 1))
# ggsave(file.path(figures_dir, "niche_vs_age_lfc_leg.pdf"), p1)
p1

p2 <- ggplot(
lfcs,
aes(x = engraf_batch, y = lfc_engraf, color = svalue_engraf < 0.1)
) +
geom_point(size = .7, alpha = 0.3) +
geom_vline(xintercept = c(-lfc_threshold , lfc_threshold)) +
geom_hline(yintercept = c(-1, 1))
# ggsave(file.path(figures_dir, "niche_vs_age_lfc_leg.pdf"), p1)
p2

p3 <- ggplot(
lfcs,
aes(x = age_batch, y = lfc_age, color = svalue_age < 0.10)
) +
geom_point(size = .7, alpha = 0.3) +
geom_vline(xintercept = c(-lfc_threshold , lfc_threshold )) +
geom_hline(yintercept = c(-1, 1))
# ggsave(file.path(figures_dir, "niche_vs_age_lfc_leg.pdf"), p1)
p3

p4 <- ggplot(
lfcs,
aes(x = niche_batch, y = lfc_niche, color = svalue_niche < 0.05)
) +
geom_point(size = .7, alpha = 0.3) +
geom_vline(xintercept = c(-lfc_threshold, lfc_threshold )) +
geom_hline(yintercept = c(-1, 1))
# ggsave(file.path(figures_dir, "niche_vs_age_lfc_leg.pdf"), p1)
p4

table(lfcs$EAN)
000 001 010 011 100 101 110 111
15204 30 59 12 46 2 6 3
ggplot(
lfcs,
aes(x = total_batch)
) +
geom_histogram(bins=60) +
scale_x_log10() +
geom_vline(xintercept = c(3))

Write to file
write_tsv(lfcs, "leg_two_batches_moderated_lfc_1_svalue_05_marked.txt")
lfcs_hits <-
lfcs %>%
dplyr::filter(lfcs$EAN != "000")
lfcs_hits
Hits by group
target_genes_leg_001 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("001")) %>%
arrange(svalue_niche)
target_genes_leg_001
figures_dir <- "./figures/hits/leg_001"
pp1 <- map(target_genes_leg_001 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)
Saving 7 x 7 in image
target_genes_leg_011 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("011")) %>%
arrange(-baseMean)
target_genes_leg_011
figures_dir <- "./figures/hits/leg_011"
pp2 <- map(target_genes_leg_011 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)
Saving 7 x 7 in image
target_genes_leg_101 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("101"))
target_genes_leg_101
figures_dir <- "./figures/hits/leg_101"
pp3 <- map(target_genes_leg_101 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)
Saving 7 x 7 in image
target_genes_leg_111 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("111"))
target_genes_leg_111
figures_dir <- "./figures/hits/leg_111"
pp4 <- map(target_genes_leg_111 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)
Saving 7 x 7 in image
target_genes_leg_100 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("100"))
target_genes_leg_100
figures_dir <- "./figures/hits/leg_100"
pp5 <- map(target_genes_leg_100 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)
Saving 7 x 7 in image
target_genes_leg_010 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("010"))
target_genes_leg_010
figures_dir <- "./figures/hits/leg_010"
pp6 <- map(target_genes_leg_010 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)
Saving 7 x 7 in image
target_genes_leg_110 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("110"))
target_genes_leg_110
figures_dir <- "./figures/hits/leg_110"
pp7 <- map(target_genes_leg_110 %>% pull(gene_id), plot_gene, data = scaled_tpm_leg_long, figures_dir)
Saving 7 x 7 in image
genes_to_plot <- c("Pax7", "Myf5", "Myod1")
figures_dir <- "./figures/additional"
map(lfcs %>%
dplyr::filter(gene_name %in% genes_to_plot) %>% pull(gene_id),
plot_gene,
data = scaled_tpm_leg_long,
figures_dir)
Saving 7 x 7 in image
[[1]]
[[2]]
[[3]]



meta_df_leg <-
as.data.frame(colData(dds_leg)[, c("age", "time", "batch")])
dt_plot_leg <-
scaled_tpm_leg %>%
mutate_if(is.numeric, ~ log1p(.))
figures_dir <- "./figures/hits"
plot_heatmap(lfcs_hits, title = "log2_scaled_tpm_leg_hits", meta_df_leg, dt_plot_leg, show_genes = F, figures_dir)

figures_dir <- "./figures/hits"
plot_heatmap(target_genes_leg_001, title = "log2_scaled_tpm_leg_hits_001", meta_df_leg, dt_plot_leg, show_genes = T, figures_dir)

figures_dir <- "./figures/hits"
plot_heatmap(target_genes_leg_100, title = "log2_scaled_tpm_leg_hits_100", meta_df_leg, dt_plot_leg, show_genes = T, figures_dir)

figures_dir <- "./figures/hits"
plot_heatmap(target_genes_leg_010, title = "log2_scaled_tpm_leg_hits_010", meta_df_leg, dt_plot_leg, show_genes = T, figures_dir)

svalue_thresh <- 0.05
figures_dir <- "./figures/volcanos"
plot_volcano(lfcs, "lfc_niche", "svalue_niche", "niche", lfc_threshold, svalue_thresh, figures_dir)
One or more p-values is 0. Converting to 10^-1 * current lowest non-zero p-value...Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

plot_volcano(lfcs, "lfc_age", "svalue_age", "age", lfc_threshold, svalue_thresh , figures_dir)
Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

plot_volcano(lfcs, "lfc_engraf", "svalue_engraf", "engraf", lfc_threshold, svalue_thresh, figures_dir)
Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] parallel stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] AnnotationHub_3.0.2 BiocFileCache_2.0.0 dbplyr_2.1.1 RColorBrewer_1.1-2 PoiClaClu_1.0.2.1 wesanderson_0.3.6 scater_1.20.1 scuttle_1.2.1
[9] SingleCellExperiment_1.14.1 ensembldb_2.16.4 AnnotationFilter_1.16.0 GenomicFeatures_1.44.2 AnnotationDbi_1.54.1 pheatmap_1.0.12 EnhancedVolcano_1.10.0 ggrepel_0.9.1
[17] rhdf5_2.36.0 apeglm_1.14.0 glmpca_0.2.0 cowplot_1.1.1 vsn_3.60.0 DESeq2_1.32.0 SummarizedExperiment_1.22.0 Biobase_2.52.0
[25] MatrixGenerics_1.4.3 matrixStats_0.61.0 GenomicRanges_1.44.0 GenomeInfoDb_1.28.4 IRanges_2.26.0 S4Vectors_0.30.2 BiocGenerics_0.38.0 tximport_1.20.0
[33] forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4 readr_2.1.2 tidyr_1.1.4 tibble_3.1.6 ggplot2_3.3.5
[41] tidyverse_1.3.1
loaded via a namespace (and not attached):
[1] utf8_1.2.2 tidyselect_1.1.1 RSQLite_2.2.9 grid_4.1.2 BiocParallel_1.26.2 munsell_0.5.0 ScaledMatrix_1.0.0
[8] ragg_1.2.1 preprocessCore_1.54.0 withr_2.4.3 colorspace_2.0-2 filelock_1.0.2 ggalt_0.4.0 knitr_1.37
[15] rstudioapi_0.13 Rttf2pt1_1.3.9 labeling_0.4.2 bbmle_1.0.24 GenomeInfoDbData_1.2.6 farver_2.1.0 bit64_4.0.5
[22] coda_0.19-4 vctrs_0.3.8 generics_0.1.1 xfun_0.29 R6_2.5.1 ggbeeswarm_0.6.0 rsvd_1.0.5
[29] locfit_1.5-9.4 bitops_1.0-7 rhdf5filters_1.4.0 cachem_1.0.6 DelayedArray_0.18.0 assertthat_0.2.1 vroom_1.5.7
[36] promises_1.2.0.1 BiocIO_1.2.0 scales_1.1.1 beeswarm_0.4.0 gtable_0.3.0 beachmat_2.8.1 ash_1.0-15
[43] affy_1.70.0 rlang_1.0.0 genefilter_1.74.1 systemfonts_1.0.3 splines_4.1.2 rtracklayer_1.52.1 extrafontdb_1.0
[50] lazyeval_0.2.2 broom_0.7.12 BiocManager_1.30.16 yaml_2.2.2 modelr_0.1.8 backports_1.4.1 httpuv_1.6.5
[57] extrafont_0.17 tools_4.1.2 affyio_1.62.0 ellipsis_0.3.2 jquerylib_0.1.4 Rcpp_1.0.8 plyr_1.8.6
[64] sparseMatrixStats_1.4.2 progress_1.2.2 zlibbioc_1.38.0 RCurl_1.98-1.5 prettyunits_1.1.1 viridis_0.6.2 haven_2.4.3
[71] fs_1.5.2 magrittr_2.0.2 reprex_2.0.1 mvtnorm_1.1-3 ProtGenerics_1.24.0 evaluate_0.14 mime_0.12
[78] hms_1.1.1 xtable_1.8-4 XML_3.99-0.8 emdbook_1.3.12 readxl_1.3.1 gridExtra_2.3 compiler_4.1.2
[85] biomaRt_2.48.3 bdsmatrix_1.3-4 maps_3.4.0 KernSmooth_2.23-20 crayon_1.4.2 htmltools_0.5.2 later_1.3.0
[92] tzdb_0.2.0 geneplotter_1.70.0 lubridate_1.8.0 DBI_1.1.2 proj4_1.0-11 MASS_7.3-55 rappdirs_0.3.3
[99] Matrix_1.4-0 cli_3.1.1 pkgconfig_2.0.3 GenomicAlignments_1.28.0 numDeriv_2016.8-1.1 xml2_1.3.3 annotate_1.70.0
[106] bslib_0.3.1 vipor_0.4.5 XVector_0.32.0 rvest_1.0.2 digest_0.6.29 Biostrings_2.60.2 rmarkdown_2.11
[113] cellranger_1.1.0 DelayedMatrixStats_1.14.3 restfulr_0.0.13 curl_4.3.2 shiny_1.7.1 Rsamtools_2.8.0 rjson_0.2.21
[120] lifecycle_1.0.1 jsonlite_1.7.3 Rhdf5lib_1.14.2 BiocNeighbors_1.10.0 viridisLite_0.4.0 limma_3.48.3 fansi_1.0.2
[127] pillar_1.6.5 lattice_0.20-45 ggrastr_1.0.1 KEGGREST_1.32.0 fastmap_1.1.0 httr_1.4.2 survival_3.2-13
[134] interactiveDisplayBase_1.30.0 glue_1.6.1 png_0.1-7 BiocVersion_3.13.1 bit_4.0.4 sass_0.4.0 stringi_1.7.6
[141] blob_1.2.2 textshaping_0.3.6 BiocSingular_1.8.1 memoise_2.0.1 irlba_2.3.5